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Abstract 

A kinetic approach for the evolution of ultracold neutral plasmas including interionic correlations 
and the treatment of ionization/excitation and recombination/deexcitation by rate equations is 
described in detail. To assess the reliability of the approximations inherent in the kinetic model, 
we have developed a hybrid molecular dynamics method. Comparison of the results reveals that 
the kinetic model describes the atomic and ionic observables of the ultracold plasma surprisingly 
well, confirming our earlier findings concerning the role of ion-ion correlations [Phys. Rev. A 68, 
010703]. In addition, the molecular dynamics approach allows one to study the relaxation of the 
ionic plasma component towards thermodynamical equilibrium. 

PACS numbers: 52.20.-j, 32.80.Pj, 52.25.Dg, 52.65.Ww 



I. INTRODUCTION 



Recent experiments have produced ultracold neutral plasmas from a small cloud of laser- 
cooled atoms confined in a magneto-optical trap . In one type of experiments 
nun , a plasma was produced by photoionizing laser-cooled Xe atoms with an initial 
ion temperature of about 10 /iK. By tuning the frequency of the ionizing laser, the initial 
electron energy E c could be varied corresponding to a temperature range IK < E e /k-Q < 
1000 K, and the subsequent expansion of the plasma into the surrounding vacuum was 
studied systematically. In a complementary type of experiment jj, 0, li| , ultracold Rb and 
Cs atoms were laser-excited into high Rydberg states rather than directly ionized. In these 
experiments, also the spontaneous evolution of the Rydberg gas into a plasma has been 
observed. The time evolution of several quantities characterizing the state of the plasma, 
such as the plasma density , the degree of ionization jj, |5|, |6( or the energy-resolved 
atomic level population [3] have been measured using various plasma diagnostic methods. 

These experiments, which have paved the way towards an unexplored regime of ionized 
gases, give rise to new phenomena in atomic physics as well as in plasma physics. Hence, a 
number of different theoretical approaches have been formulated to cover different aspects 
of these experiments B, 3 3, H U, 12, 3- 



An important issue is the question whether the plasma produced would be strongly cou- 
pled or not. The correlation strength is determined by the Coulomb coupling parameter 
r = e 2 /{ak-oT) with the Wigner-Seitz radius a A plasma is called "strongly coupled" if 
r>l, i.e. if the Coulomb interaction between the plasma particles largely exceeds the ther- 
mal kinetic ener gy. In this case, interesting ordering effects such as Coulomb crystallization 
can be observed [jjj, For the initial conditions of the NIST experiments Q,Q,0|, however, 
the development of equilibrium electron-electron correlations leads to a rapid heating of the 
electron gas, which prevents the electron Coulomb coupling parameter T e from exceeding 
unity 0. The same has been argued for an ion plasma in Q|. Since the electron dynamics 
proceeds on a much smaller timescale than the ion motion, in [3, [| the electron heating 
could be studied for the early stage of the plasma evolution only, where the ionic component 
does not show dynamical effects. On the other hand, ion heating has only been studied in 
;he framework of a model system, consisting of a homogeneous gas of Debye-screened ions 
~\\ . such that the influence of the subsequent expansion could not be explored. 
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The first quantitative comparison with experimental observations has been given in [lo| . 
with the plasma dynamics modeled within a hydrodynamical approach and ionization, exci- 
tation and recombination treated by a separate set of rate equations. Since this model does 
only account for the mean-field potential created by the charges, it cannot describe effects 
of particle correlations. However, it has also been shown there that the electronic Coulomb 
coupling parameter T e does not exceed a value of ~ 0.2 during the plasma expansion due to 
heating by three-body recombination. Thus, the influence of electron-electron correlations 
on the dynamics of the plasma could be neglected. On the other hand, three-body recombi- 
nation does not influence the ionic temperature, so that the ions can heat up only through 
correlation heating (and energy exchange with the electrons, which, however, is very slow). 
Since the ionic temperature was set to zero in |10], the role of ion-ion correlations could not 

n 

be explored. In a preliminary study [13J , we showed that they indeed change the evolution of 
the system quantitatively, though not qualitatively. In the following, we will give a detailed 
account of the kinetic model used in |l3j and of all relevant ingredients. We will also develop 
a hybrid molecular dynamics (H-MD) approach which treats the electronic plasma compo- 
nent in an adiabatic approximation while the ions are fully accounted for. Such an approach 
permits the description of situations where the ions are strongly coupled [tgI . Il^ , which is 
clearly beyond the capabilities of the simple kinetic model. Nevertheless, for the typical 
situations realized in the experiments flflQ, comparison of the two theoretica. approaches 
yields very good agreement, corroborating our findings reported earlier (13J and establishing 
firmly that one can capture the relevant physics with the relatively simple kinetic approach. 



II. THEORETICAL APPROACH 

Our kinetic approach is similar to the one of Q|. The main difference is the inclu- 
sion of ion-ion correlations (IC) which will be described in detail below. Briefly, a set of 
kinetic equations is formulated for the evolution of the plasma fsubsection lll A|) . while ioniza- 
tion/excitation and recombination/deexcitation are taken into account on the basis of rate 
equations ( subsection IIIB|) . In order to test the applicability and accuracy of this model, 
we have developed a less approximative and more flexible but computationally much more 
demanding approach. It uses molecular dynamics for the ionic motion while the electron 
component is treated as a fluid assuming a quasi-steady state (subsection III Cj) . 
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A. Kinetic description 



Starting from the first equation of the BBGKY hierarchy, the evolution equation for the 
one-particle distribution function f a (r,v,t) of the free plasma charges is obtained as 

m a (d t + v<9 r ) f a (r, v,t) = ^2 / [drVapi*, r')] d v f af3 (r, v, r', v', t) dr'dV , (1) 

where a, (3 label the particle species (e,i for electrons and ions, respectively), fap(r, v, r', v', t) 
denotes the two-particle distribution function for the corresponding particle species and 
fa/3 = la^p/lf — r '\ is the Coulomb interaction potential between the charges q a and qp. 
Electron-electron correlations are very small during the plasma expansion, since the electrons 
will quickly heat up due to three-body recombination and the additional heating due to 
correlation effects is small in comparison [10]. Hence, we neglect electron-electron as well 
as electron-ion correlations, leaving only IC as a possible influence on the plasma dynamics 
beyond the mean-field level. On this level of approximation the ion kinetic equation can be 
written as 

m (d t + vd r - ^-dv^j fi = J (<9 r y?ii) 9 v w u (r, v, r', v') dr'dv' , (2) 
where the function 

I0ii(r, v, r, v') = /ii(r, v, r, v) - /i(r, v)/i(r', v') (3) 

contains the contributions of IC to the two-particle distribution function and (p is the mean- 
field potential created by all plasma charges. Since m e jm- x 1, the relaxation timescale 
of the electrons is much smaller than the timescale of the plasma expansion under typical 
experimental conditions |lj. Thus, we may safely apply an adiabatic approximation for the 
electron distribution function, assuming a local Maxwellian distribution 

/e(r,v)ocexp(|g) exp (-^r) , (4) 

where T e is the electron temperature. Eq. (j3J) together with a quasineutral approximation 
[rol l allows one to express the mean-field potential in terms of the ionic density P\ = j f\dv, 
resulting in 

d r ip = k B T e ^- . (5) 

Pi 
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Using Eq. (jSJ), the following evolution equations for the second moments of the ion distribu- 
tion function are derived from Eq. (J2J) 

d t (r 2 ) = 2(rv) (6a) 

Y dt < rv > = Y ( v ") + l k * T * + \ N i l j Pi(r)rF ii (r) dr (6b) 

^±d t (v 2 ) = Nr 1 (^k B T c J Pi d T u dr- J v(9 r ^ ii )w ii (r, v, r', v') drdvdr'dv') (6c) 
where (r 2 ) = N^ 1 f r 2 /i(r, v)drdv etc. The "correlation force" F^r) is given by 

Fh (r) = - J few) Pi (r') 9ii (r, r') dr' , (7) 

where the spatial correlation function ga is defined by pi(r)pi(r')<?ii(r, r') = 
J Wii (r, v, r', v') d\dV and u(r) = j vfidw is the hydrodynamical drift velocity of the 
plasma. With the help of the second kinetic ion equation of the BBGKY hierarchy, the 
last term on the right-hand side of Eq. (JoT^) can be written as 

— f v(9 r <y9 ii )w i i drdvdr'dv' = —d t f faWa drdvdr'dv' 

= -d t Un , (8) 

where 

Uii = wJ V?nPi(r)Pi( r ')2ii(r, r') drdr' (9) 
is the average correlation energy per ion. Hence, Eq. (j6^|) reflects energy conservation for 
the ion subsystem. The evolution of the hydrodynamical velocity u is determined by 

m iPi [d t u + (u • <9 r ) u] = -k B T c d r pi - d r P th;i - p^n (10) 

where P t h,i = ^ / ( v — u ) 2 /i drdv is the thermal ion pressure. As shown in the appendix, 
in the framework of a local density approximation, i.e. by assuming that g$ only depends 
on the distance |r — r'| and on the densities at the two coordinates r and r', and that the 
ionic density p^ varies slowly on the lengthscale where g is significantly different from zero, 
the total correlation energy can be approximated by the well-known LDA expression 

Uu = Nr 1 J UiiPi dr (11) 

while the correlation force is found to be 

F„ = -i^ + ^W. • (12) 
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where «u(r) is the correlation energy of a homogeneous plasma of density Pi(r), 

«H(r) = yp»(r) | ^-^<lx . (13) 

If IC are neglected in Eq. (J2J), the kinetic equation exhibits the following selfsimilar solution 

f r 2 \ ( rrn (v - u) 2 
/i oc exp -— — exp 



2cr 2 y M 2fc B Ti 
u = 7r (14) 

which corresponds to the initial state of the experiments under consideration. As soon as IC 
are taken into account via the correlation pressure piF;; in Eq. (JTUj), however, Eqs. (fTlj) are 
no longer exact solutions of Eq. (J2J). Using Eq. (|12p. the last term on the right-hand side of 
Eq. (|T0j) can be rewritten as piFa = — \ 9 ^ lpl ' > d r p\. Interpre ting this term as a local nonideal 



by averaging the differential 
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pressure, an equation for the parameter 7 was derived in 
equation for |u|/r obtained by inserting Eq. (jTj| into Eq. (fTUJ) over the plasma volume. 
Obviously, this treatment is not unique. Since, as discussed above, the ansatz (JTljl does not 
solve Eq. (fTU|) exactly, multiplying Eq. (fTU|) by different functions of r and averaging over 
the plasma volume will lead to slightly different evolution equations for the parameter 7. 
Here, supported by a comparison with our numerical results from the MD simulations, to 
be discussed below, we assume as in [.13] that the functional form of the hydrodynamical 
quantities of Eqs. (fTf|) is not altered by the inclusion of IC, while the dynamics of the 
parameters appearing in Eqs. (JUjl is determined from the equations © for the moments 
of the distribution function. Clearly, such an approximation can not be a priori justified. 
Hence, it must be validated a posteriori by comparison with more sophisticated methods 
which do not rely on a reduction of the plasma description to a few macroscopic parameters. 

With this procedure, we arrive at the following set of equations for the width a of the 
plasma cloud, its expansion velocity 7r as well as ionic and electronic temperature 7] and 
T • 



d t a 2 = 2 7 a 2 (15a) 



^2 



= n „ ■„ , 3 ^ _ y, (15b) 

d t k B Ti = -27M - pU a - p t U n (15c) 
d t k B T e = -2jk B T e . (15d) 
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The last equation (jl5d|) has been derived from the electron kinetic equation by making use 
of the quasineutrality condition. The set of Eqs. (fT5|) slightly differs from that presented in 



Q 



13] where W c = -j^ J Pi 9 ^^ dr had been used instead of Ua in Eq. (j!5b|) . A comparison 
with our numerical MD results shows that Eqs. 1)1 5|) yield a slightly better quantitative 
agreement, while the principal influence of IC on the plasma dynamics, which has been 
partly discussed in 13], is the same. Eqs. (|T5j) provide a transparent physical picture of 
the expansion dynamics. First, Eq. (jl5dj) together with (j!5a|) reflects the adiabatic cooling 
of the electron gas, i.e. T e a 2 = const. The ion temperature, on the other hand, is not 
only affected by the adiabatic cooling, expressed by the first term in Eq. 1)1 5c|) . but also 
changes due to the development of IC, which is taken into account by the last term in Eq. 
(|15cj) . Furthermore, these correlations reduce the ion- ion interaction and therefore lead to 
an effective negative acceleration, expressed by the Un/3-term in Eq. (jl5b|) . in addition to 
the ideal thermal pressure. This contribution, which corresponds to the average nonideal 
pressure known from homogeneous systems llfij . also leads to an effective potential in 
which the ions move. As they expand in this potential, the thermal energy changes due to 
energy conservation, as expressed by the second term on the right-hand side of Eq. I)15cj) . 
Finally, combining eqs. ()15|) yields a second integral of motion, namely the total energy of 
the plasma 

E tot = - {k B T e + fc B Ti) + -m i7 V + U a . (16) 

Although the set of equations ()15|) determines the time evolution of all relevant macroscopic 
plasma parameters, namely its width, expansion velocity, electron and ion temperature, it 
is not a closed set since an evolution equation for the correlation energy Ua which enters 
Eq. ()15c|) is missing. Initially, the plasma is completely uncorrelated, so that Ua(t = 0) = 0. 
However, the initial state corresponds to a non-equilibrium situation, and the plasma will 
relax towards thermodynamic equilibrium, thereby building up correlations. A precise de- 
scription of this relaxation process in the framework of a kinetic theory is rather complicated 
and requires a considerable numerical effort j^J. We therefore employ a linear approxima- 
tion for the relaxation of the two-particle correlation function, the so-called correlation-time 



approximation 



2l|, 



dw n (r,v,r',v';t) _ Wq(r, v, r', v'; t) - w-?(r, r'; t) ^ 



dt r< 



corr 



Here, r c , 
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is the characteristic timescale for the relaxation of particle correlations and 



u^ q is the equilibrium pair correlation function, which in our case still depends on time via the 
evolving one-particle distribution function since the plasma is freely expanding. As shown 
in [23! the correlation time r corr can be well estimated by the inverse ionic plasma frequency. 



Hence, in our calculations we set r corr = u~\ = a/ m-J (4ice 2 p i ), where pi = N-J (Ana 2 ) 3 / 2 is 
the average ionic density of the plasma. Such a linear approximation is good only for small 
deviations of from its equilibrium form. Clearly, this is not the case in the initial stage 
of the gas evolution. However, after the initial phase of correlation heating the system stays 
very close to its slowly changing local equilibrium, and one may expect Eq. (JT7j) to yield 
good results. Under the same conditions that lead to Eqs. (fTTJl and (JI2J), one easily verifies 
that Eq. flT7D leads to 

d t Uz « - " U » , (18) 

''"corr 

where = N-' 1 J piU^dr and u eq (r) is the correlation energy per particle of a homogeneous 
one-component plasma in local equilibrium. This quantity has been studied intensively in 
the past, and approximate analytical formulae are available in the literature Here, 
we adopt the interpolation formula from 2^] 

<" = teT ' r3/2 faftr + itr) ■ (19) 

with A x = -0.9052, A 2 = 0.6322 and A 3 = -y/3/2 - A x j ' y/M, which yields an accurate 
interpolation between the low-r Abe limit and the high-r behavior obtained by Monte Carlo 
and MD simulations. It should be noted that in the present situation uu depends on time 
since the plasma expands. Hence, T and with it the thermodynamical equilibrium change 
in time. 

The set of equations (|15|) describes the evolution of the plasma part of the system, i.e. a 
system of Ni ions and electrons. Due to ionization and recombination events occurring during 
the plasma expansion (discussed in detail in the following subsection), this number N u and 
hence also the total mass M = N i m 3jtom , is not constant over the course of the evolution. 
However, such a treatment completely neglects the influence of the bound Rydberg atoms 
on the dynamics. One may argue that they do not influence the plasma evolution since 
they do not interact with the ions or electrons by Coulomb interaction. On the other 
hand, a Rydberg atom may carry a significant amount of kinetic energy, gained from the 
acceleration by the electron pressure before its formation by three-body recombination. In 
a simple approximation, we assume equal hydrodynamical velocities and density profiles for 
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the ions and atoms, in order to account for this effect. This implies that the expansion of 
the neutral Rydberg atoms can be taken into account by replacing the mass N\m\ of the 
ions by the mass of the total system [N\ + iV a )mi, where N a is the number of atoms. We 
therefore replace the ion mass m,\ by an effective mass (1 + NjN^rrii in Eq. (fl5b"jl . The 
quality of this approximation can, of course, also be checked by comparison with the H-MD 
description, see below. 



B. Ionization and Recombination 

As demonstrated in , a satisfactory description of the dynamics of an ultracold plasma 
can be achieved by combining a hydrodynamic treatment of the plasma evolution with 
rate equations accounting for inelastic collisions between the plasma particles and Rydberg 
atoms. The rate equation for the change of density of Rydberg atoms in a state with 
principal quantum number n reads 



^2 [ K {v, n)p a (p) - K(n, p)p a (n)] + p c [R(n)p c pi - 7(n)p a (n)] , (20) 



where K (p, n) is the rate coefficient for electron impact (de)excitation from level p to level n, 
and R(n) and I(n) describe three-body recombination into and electron-impact ionization 
from level n, respectively. The rate coefficients K, R and / have been taken from the classic 
wo^ofMa^andKeckQ. AddHioaa! processes, such as, e.g., ionization by black- 
body radiation or from dipolar atom-atom interactions, are easily included in Eq. ()20)1 if the 
corresponding rates are available. Such processes are essential for a description of the early 
stages of the evolution of a system starting with a Rydberg gas Q, 0, 0], but are of minor 
importance in situations starting from a pure plasma. 

In this framework, the evolution of the system is obtained by solving Eqs. ()15aJ) - ()15c|) 
together with Eq. ()20|) while the electron temperature is now obtained from the modified en- 
ergy conservation relation E tot +E a = const, instead of Eq. (|15dj) . where E a = —1Z J2 n N^n' 2 
is the total energy of the Rydberg atoms and TZ = 13.6 eV. 
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C. Hybrid molecular dynamics treatment 

As we will show in section III, the kinetic description of the previous subsections is 
able to describe the plasma dynamics to a surprisingly large extent. However, one of the 
main motivations of this work is the study of the role of IC, which are incorporated in 
the model only in an approximate way. To assess their influence on the dynamics reliably, 
a more sophisticated approach is required, e.g. molecular dynamics simulations which fully 
incorporate the ionic interactions. However, a full MD simulation of both, electrons and ions, 
is computationally very demanding, and only the very early stage of the system evolution 
can be described in this way |?| • On the other hand, as argued above, electronic correlations 
are not important for the plasma dynamics, so that only IC have to be accounted for in 
full while the influence of the electrons on the dynamics may be treated on a mean-field 
level. Moreover, we have seen that the timescale of equilibration of the electronic subsystem 
is orders of magnitude shorter than that of the ionic subsystem and the timescale of the 
plasma expansion. This observation led us to use an adiabatic approximation in subsection 
111 Al where the electrons are assumed to equilibrate instantaneously, assuming a Maxwellian 
velocity distribution with a well-defined temperature and a spatial profile determined from 
the total mean-field potential of the plasma charges. The clear separation of timescales 
suggests that this adiabatic approximation is well justified, hence we will keep it in the 
following. Consequently, we have developed a hybrid approach where the electrons are 
treated on a hydrodynamical level as in the kinetic description above, while the ions are 
propagated individually with their mutual interaction and the influence of the electrons on 
the ions enters via the electronic mean-field potential. This hybrid approach permits the 
use of much larger timesteps in the propagation of the system, since the electronic dynamics 
needs not to be followed in detail but only the ionic motion has to be resolved in time. 
Consequently, the evolution of the system can be followed over the experimental timescales. 
Furthermore, the approximate treatment of IC in the kinetic model of subsection III Al can 
be tested. Finally, beyond the scope of the present work, we have shown that the 
present H-MD approach can describe situations where the ionic plasma component is so 
strongly coupled that crystallization of the ions sets in. Such a scenario is clearly beyond 
the capabilities of a kinetic approach. 

As discussed above, the electrons are still treated as a fluid, while we lift the quasineutral 
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approximation by calculating the resulting mean-field potential from the Poisson equation 

Atp = Aire 2 (p e — p\) . (21) 



However, using Eq. (jlj) poses a conceptual difficulty |ll| since the mean-field potential ap- 
proaches a finite value at large distances and therefore leads to a non-normalizable electron 
density. This problem, which has been discussed for a long time in an astrophysical context 
2(| , reflects the fact that a substantial fraction of the electrons indeed escapes the finite po- 
tential barrier at long times during the relaxation process until the total kinetic energy of all 
electrons is less than the height of the potential well. On the timescales under consideration, 
however, typically only a small amount of the electrons escapes the plasma volume, until 
the resulting charge imbalance becomes large enough to trap the remaining electrons, which 
quickly reach a quasi-steady state forming a temporarily quasineutral plasma in the central 
region. We account for this electron loss by determining the fraction of trapped electrons 
from the results of Ref. 

The corresponding steady-state distribution, derived for the study of globular clusters, is 
of the form |27] 

p e (x exp {j^Tj^J J exp (- x ) x3/2 dx i ( 22 ) 
where x — m e^esc/(2^B^c) with the velocity v esc (r) necessary to escape from a given position 
in the plasma. In the present case, the potential can have a non-monotonous radial space 
dependence and the escape velocity has to be defined as 

TTl 

"TT^esc ( r ) = max [<f (r') - <p (r)], (23) 

Z r'>r 

in contrast to astrophysical problems where one only has a single sign of "charge" and 
m e v 2 sc /2 = —<p |27[. For a given electron temperature T e and ion density p\ the electron 
density is found by numerical iteration of Eqs. (}2*Tj) . (j2^j) and (|2^j) until selfconsistency is 
reached. 

Knowledge of the electron density then permits propagation of the ions in the electron 
mean-field Aip c = 47re 2 p e and the full interaction potential of the remaining ions, 

m^ = d r ^ c + e 2 J2 , rj ~ r ,3 • (24) 
k l r i ~ r k\ 

The numerical solution of the ion equations of motion represents the most time consuming 
part of the plasma propagation. In general, for N propagated particles, the corresponding 
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numerical effort scales with N 2 rendering a treatment of large particle numbers difficult. 
In order to simulate particle numbers relevant to the experiments, we have adapted a hi- 
erarchical treecode originally designed for astrophysical problems, first described in 
This method provides a numerically exact solution of the ion equations of motion Eq. (|24j) . 
while the numerical effort grows only as N In N with increasing N. More details about the 
numerical procedure can be found, e.g., in 2^ |. 

In the framework of the kinetic model introduced in section III Al the influence of IC 
on the system evolution can be singled out by comparison with the solution of the corre- 
sponding equations with Ua = 0. In order to make an analogous comparison also for the 
MD simulations, we have performed calculations propagating the ions in the mean-field po- 
tential created by all charges. Technically, the mean-field potential is represented using a 
test-particle method, widely used for various problems in plasma physics (see, e.g., j^])- 



III. RESULTS AND DISCUSSION 



We will discuss the evolution of a plasma initially consisting of N e = 37500 electrons 
and iV; = 40000 ions with an average density of 10 9 cm -3 at a rather low electronic kinetic 
energy E e = ?>k&T e /2 = 20 K, comparing the results from the kinetic model and our MD 
simulation. Thereby, we put special emphasis on the role of IC. 



A. Global aspects of plasma expansion and recombination 

The general macroscopic behavior of the system has been described before in several 
publications, experimentally as well as theoretically , [ill [l]] . The plasma cloud slowly 
expands due to the thermal pressure of the electrons, leading to adiabatic cooling of the 
electrons as well as partial recombination into bound states (figures ^ and |2J) • The amount 
of recombination and its influence strongly depends on the initial electron temperature and 
density. If the electrons are too hot (about 50 K for typical experimental densities of 10 9 
cm" 3 ), recombination is strongly suppressed and the system dynamics is well described by 
the results of |3] obtained for the collisionless plasma expansion 
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FIG. 1: Electronic temperature T e (t) for an expanding plasma of 40000 ions with an initial 
average density of 10 9 cm -3 and an initial electron kinetic energy of 20 K, obtained from the H-MD 
simulation (a) and the kinetic model (b), with (solid) and without (dotted) the inclusion of IC. 
The inset shows the ratio of the electron temperatures obtained from the H-MD simulation and 
the kinetic model. 

1. Temporal evolution of the electronic temperature 

For the lower electron temperatures considered here, as can be seen in Fig. there is an 
initial increase of the electron temperature due to electron heating by three-body recombi- 
nation and subsequent deexcitation of the formed Rydberg atoms. At low initial electron 
energies this heating drastically increases the electron temperature and thus accelerates the 
plasma expansion jjLjj], which explains the enhanced expansion velocity observed in In 
contrast to this recombination heating of the electrons, the inclusion of IC only slightly 
changes the expansion dynamics, as seen in Fig. Eby comparing the solid and dotted lines. 
As shown in the inset of Fig. HJi, the electron temperature obtained from the H-MD sim- 
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ulation and the kinetic model differ by at most 8% during the first few microseconds of 
the plasma expansion, while the agreement becomes even better at later times. Moreover, 
the faster decrease of the electron temperature due to the inclusion of IC, predicted by the 
particle simulations, is quantitatively reproduced by the much simpler kinetic model. 

Hence, the simple evolution equations (|15|) are sufficient to clarify the role of IC in 
the expansion dynamics. According to Eq. ()15c|) . the development of IC quickly heats 
up the plasma ions to roughly —§£4 since the expansion of the plasma is still negligible 
during this initial stage. Thereby, the negative correlation energy term ^Ua in Eq. (jl5b|) 
is overcompensated, leading to a faster expansion of the plasma. As a consequence of the 
quicker expansion, the electron temperature decreases somewhat faster than without the 
inclusion of IC. With Eq. (j!5b|) . the importance of this effect can be estimated by comparing 
the thermal electron energy /cb^c to the net ion contribution — ^Ua in the numerator of the 
first term on the right-hand side of Eq. (|15b|) . Estimating the correlation energy by e 2 /a, 
it follows that the total pressure driving the plasma expansion is enhanced by a factor of 
roughly 1 + r e /3, which only slightly changes the expansion dynamics since the electrons 
are known to be weakly coupled over the whole observation time |lfl| . 



2. Formation of Rydberg atoms in time 

The number of recombined atoms is influenced more strongly by IC (Fig. EJ). During the 

evolution of the system, Rydberg atoms are constantly formed by three-body recombination 

and re-ionized by the free electrons in the plasma. As shown in Fig. |2t, for the current set 

of parameters about 7000 Rydberg atoms are present in the system after 40 /is, while the 

kinetic model yields about 7900 atoms at the same instant of time (Fig. |2b). This number 

is small compared to the size of the whole system, nevertheless it is large enough that the 

recombined atoms can be detected in an experiment, and corresponding curves have indeed 

fl 

been obtained experimentally [3J|. Due to the strong temperature dependence of the total 
three-body recombination rate, which is proportional to T e -9 ^ 2 (^J, the slight decrease of 
the electron temperature due to the faster expansion, caused by the correlation heating of 
the ions, considerably affects the recombination behavior of the plasma. While there is 
an overall shift between the atom number obtained from the particle simulations and the 
kinetic model, both the kinetic model and the H-MD simulation yield an increase of the 
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FIG. 2: Number N a (t) of recombined atoms obtained from the H-MD simulation (a) and the 
kinetic model (b) for the same parameters as in Fig. ^ The solid line shows the result taking into 
account the IC while the dotted line is obtained from the mean-field treatment. 

atom number of about 10% at t — 40 fis (Fig. EJ), compared to a mean-field treatment of the 
ion dynamics. Thus, the H-MD simulation corroborates our previous findings 

Additional insight into the recombination process can be gained from a closer look at the 
distribution of bound Rydberg states. Figure El shows the population of levels with principal 
quantum number n for three different times, corresponding to different stages of the plasma 
expansion. Initially, Rydberg states of moderate excitation are populated, due to a relatively 
high electron temperature (Fig.Et). At later times, higher excited bound states are formed in 
the course of the plasma expansion (Figs.Eb andEt), since the maximum principal quantum 
number for recombination n max = \flZ/ (2kBT e ) [25j increases as the electron temperature 
drops down. Moreover, the deeply bound states formed at earlier times are also not subject 
to electron-impact excitation and deexcitation anymore since the thermal velocity of the 
impacting electrons has become too small. Thus, as becomes apparent by comparing Fig. 
Efc> with Eh, the deeply bound states (n ^5 30) remain basically untouched, while higher and 
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FIG. 3: Population of bound Rydberg states with principal quantum number re, after t = 1.5 /is 
(a), t = 6 /is (b) and t = 40 /xs (c). The vertical bars represent the H-MD calculation, the solid 
line the kinetic model. The dashed curve in (c) shows the kinetic model neglecting IC. Initial-state 
parameters are the same as in Fig. ^ 



higher states "freeze out" as the plasma expands. As may be anticipated from Fig. El IC 
mainly affect the later stages of the plasma evolution. Hence, the inclusion of IC alters the 
population of these higher lying states, as shown in Fig. Et- Since these states have small 
binding energy, they contribute little to the total kinetic energy of the plasma subsystem. 
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FIG. 4: Spatial densities p\ (solid) and p a (dashed) of the ions and atoms, respectively, at t = 3 /xs 
(a) and t = 31.3 fis (b), compared to the Gaussian profile assumed for the kinetic model (dotted). 
Additionally, p\ obtained from the particle simulation using the mean-field interaction only is shown 
as the dot-dashed line in (a). Initial-state parameters are the same as in Fig. ^ 

This is the reason why the effect of IC is visible in the distribution of Rydberg states, but 
not in the macroscopic expansion dynamics of the plasma, reflected, e.g., by the asymptotic 
expansion velocity measured in Q]. 



B. Spatially resolved plasma expansion and relaxation 

While the time evolution of global, i.e. space-averaged, observables of the plasma is very 
well described by the kinetic model, one may expect discrepancies compared to the MD 
simulations when looking into the spatially resolved plasma dynamics. We will assess these 
discrepancies quantitatively in the following. 



17 



-3QL 



100 200 300 400 500 

r [jam] 




1000 2000 3000 4000 5000 
r [jam] 



FIG. 5: Hydrodynamic velocity u(r) of ions (full circles) and atoms (open circles) at the same 
times as in Fig. 0J compared to the straight-line assumption of the kinetic model. The dashed 
line in (a) shows the result of the particle simulation using a mean-field treatment of the ion-ion 
interaction. 



1. Evolution of the particle densities 

In the derivation of the kinetic equations (j!5|) . we have assumed that the analytical form of 
the ionic density p x remains invariant during the evolution of the system and, moreover, that 
the atoms will have the same distribution. As the plasma expands, the spatial profile of the 
ions must deviate from its original Gaussian shape 11]. This is mainly due to deviations from 
quasineutrality, e.g. deviations from the linear space dependence of the outward directed 
acceleration, at the plasma edge. The influence of the nonlinear correlation pressure on 
the density profile is of minor importance, as can be seen by comparing the solid and dot- 
dashed line of Fig. in the inner plasma region. As known from earlier studies of expanding 
plasmas, based on a mean-field treatment of the particle interactions a sharp 

spike develops at the plasma edge, shown by the dot-dashed line in Fig. 0^,. At later times, 
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this spike decays again when the maximum of the hydro dynamic ion velocity passes the 
position of the density peak, so that the region of the peak is depleted. Ultimately, at 
long times, the plasma approaches a quasineutral selfsimilar expansion |32j ]. From Fig. 0^ 
it becomes apparent that with IC the peak structure is less pronounced than in mean-field 
approximation. This is due to dissipation caused by ion-ion collisions which are fully taken 
into account in the H-MD simulation. As shown in 32], by adding an ion viscosity term 
to the hydrodynamic equations of motion, dissipation tends to stabilize the ion density 
and prevents the occurrence of wavebreaking which was found to be responsible for the 
diverging ion density at the plasma edge in the case of a dissipationless plasma expansion. 
Furthermore, the initial correlation heating of the ions largely increases the thermal ion 
velocities leading to a broadening of the peak structure compared to the zero-temperature 
case. 

Apart from the deviations at the plasma edge, the ionic density is rather well reproduced 
by the Gaussian approximation for the spatial distribution. In particular, there is good 
agreement between the rms-radii obtained from the MD simulation and the kinetic model. 
On the other hand, the spatial distribution of atoms significantly deviates from that of the 
ions even at relatively early times due to the nonlinear density dependence of the collision 
rates in Eq. (|20|). However, as also stated in jn|, the total number of atoms is too small to 
significantly influence the macroscopic expansion of the system. 



2. Spatial dependence of the radial velocities 

Another assumption used in the derivation of the kinetic model is the proportionality 
of the hydrodynamical expansion velocity to the distance from the center of the plasma 
cloud, u = 7r, both for the ions and the atoms. In order to check this assumption, we have 
calculated the radial velocity component v r = vr/r of each particle, which is plotted as a 
function of the radial distance from the plasma center in Fig. |SJ At an early time the velocity 
distribution is spread out about its mean value predicted from the kinetic approach due to 
the finite ionic temperature. Note that the expansion is slower near the plasma edge due to 
the deviation from quasineutrality as discussed above. Consequently, the inner part of the 
plasma which expands more quickly will catch up with the outer rim, leading to the formation 
of the density spike seen in Fig. 0^,. In the case of the H-MD simulation the velocity spread, 
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FIG. 6: Distribution of thermal ionic velocities at t = 0.3 /xs sampled from three different regions 
of the plasma: r < 1.3a (a), 1.3a < r < 2a (b), and r > 2a (c), and from the total plasma 
volume (d). The solid lines show a fit to a Maxwell-Boltzmann distribution corresponding to the 
temperatures specified in the figure. Initial-state parameters are the same as in Fig. ^ 

caused by the initial ion heating, is of the same order of magnitude as the hydrodynamical 
expansion velocity itself, leading to a significant broadening of the density spike as discussed 
above. At later stages of the system evolution, the ions cool adiabatically due to the plasma 
expansion, and the width of the velocity distribution decreases significantly. Moreover, as 
discussed in connection with the decay of the ion density peak in Fig. Eb, the decrease of 
the ion velocities near the plasma edge apparent at early times has disappeared. 

A comparison with the result of the kinetic model equations (JT5j) shows once more that the 
H-MD simulation not only reproduces the linear radial dependence of the hydrodynamical 
velocity, but also yields a quantitative agreement between both methods. 



3. Spatial dependence of the thermal velocities 

Due to its marginal influence on the plasma expansion dynamics, the role of the ionic 
temperature Tj for the state of the system has not been addressed before. In the cold fluid 
model of lfj, [llj, 7] has been set to zero in order to follow the long-time plasma dynamics. 
However, as stated in the introduction, one of the motivations of the current type of exper- 
iments was the creation of a strongly coupled plasma. In this context, knowledge of 7] is 
essential since it directly enters the Coulomb coupling parameter which determines the state 
of the plasma. Moreover, the ionic temperature gives important insight into the relaxation 
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FIG. 7: Same as Fig. El but for t = 1.5 fj,s. 



dynamics of the plasma. For comparing the kinetic model with the H-MD calculations, the 
very definition of T; for the MD simulation requires some discussion. As discussed in section 
II, we assume a Gaussian velocity distribution, i.e. a well-defined temperature TJ, for the 
plasma ions in our kinetic model. This, of course, is an approximation since the plasma is 
not created in an equilibrium state. The total kinetic energy of the ions is a sum of the 
hydrodynamical expansion energy and a contribution due to the thermal motion of the ions. 
Since the hydrodynamical velocity is directed radially (Eq. (^3)), we determine the thermal 
energy of the ions from the average of the velocity component perpendicular to the radial 
direction 



Clearly, such an assignment of a temperature to the average velocity is only well defined if 
the ion velocities v± are distributed according to a Maxwell distribution. In order to check 
the validity of this requirement, we have sampled the ion velocity distribution v± from three 
different regions in the plasma: r < 1.3a, 1.3a < r < 2a, and r > 2a, which have been 
chosen so that each region is occupied by approximately the same number of ions. The 
resulting distributions are plotted at two different times t = 0.3 /is = 1-30;"^ and t = 1.5 fxs 
= in Figs. El and respectively. Additionally, we have fitted a Maxwell-Boltzmann 

distribution to the numerical results, formally defining a temperature in the corresponding 
plasma region. As can be seen in Fig. even at the very early stage of the plasma evolution 
the numerical data is well fitted by an equilibrium distribution in the inner plasma region, 
while there are deviations in the outer region of the plasma since the relaxation time is longer 




(25) 
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FIG. 8: Average thermal ionic energy as a function of the distance from the plasma center at 
four different times: t = 0.3 fis (solid), t = 1.5 /j,s (dashed), t = 6.0 /is (dotted) and t = 25.3 /is 
(dot-dashed) . Initial-state parameters are the same as in Fig. ^ 

due to the lower density far away from the plasma center. However, already after a relatively 
short time of t — 1.5 /is the velocity distributions are well fitted by a Maxwell-Boltzmann 
distribution in all three plasma regions (Fig. |7j). Hence, the ion thermal energy can be 
represented by a local temperature 71 (r), decreasing with growing distance from the plasma 
center as can be seen from Figs. El and UJ This is due to the fact that the initial heating 
arises from a compensation of the negative correlation energy, which is larger in the central 
plasma region where the density is higher. However, as becomes apparent by comparing 
Figs. El and [71 the thermal energy equilibrates over the whole plasma volume rather quickly 
as the system evolves. While the temperatures denned in the inner and outermost region 
deviate by a factor of eight at t = 0.3 /is, they differ by a factor of two only 1.2 /is later. 

Figure |H1 gives a more detailed account of this equilibration process. Here, the local ionic 
temperature is plotted as a function of the radial distance from the plasma center at four 
different times, where 71 (r) has been defined from the velocity average of a shell of 2000 
ions with a central shell radius r. The temperature decrease with increasing distance from 
the center as discussed above is clearly visible. Nevertheless, the ion temperature is seen 
to equilibrate rather quickly, so that the approximation of a homogeneous ion temperature, 
used in the derivation of the kinetic model in section III A| becomes better and better at later 
times. Moreover, the numerically calculated distribution of thermal velocities sampled over 
the whole plasma volume is well represented by a Maxwell-Boltzmann distribution with some 
average temperature intermediate between the temperatures of the inner and outer region, 
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respectively (Fig. 01). This shows that the Gaussian phase-space distribution assumed for 
the ions in section III Al agrees very well with the results of the MD simulation averaged over 
the spatial coordinates, even if the temperature still shows substantial inhomogeneities. 



C. Spatially averaged ionic observables 

As we have demonstrated in section UlI Al the kinetic model describes the global temporal 
evolution of the plasma including recombination quite accurately. From the detailed anal- 
ysis of the spatially resolved plasma dynamics in the previous subsection we may expect 
that the kinetic model describes spatially averaged observables, such as the kinetic energy 
of the expansion, the thermal energy, and the correlation energy of the plasma quite well. 
This is indeed the case over almost the entire evolution time as Fig. [§] demonstrates for 
the correlation energy and the thermal ion energy. Only at an early stage of the plasma 
evolution, differences between MD simulation and kinetic model are visible, showing that 
the correlation-time approximation Eq. ()17j) does not accurately describe this early phase 
of equilibration starting from a completely uncorrelated state in all details. Since the initial 
state is very far from equilibrium, the initial relaxation process is not exponential, as as- 
sumed in the correlation-time approximation Eq. ([17)1 . Rather, it is connected with transient 
oscillations of the temperature (inset of Fig. El which have been found both theoretically 

nnn n 

[la |3Ja, |2J] and experimentally [351 ] . However, the timescale of the initial ion heating as well 
as the maximum temperature are well reproduced by the simple model. After the system has 
come sufficiently close to local equilibrium, the quality of the correlation-time approxima- 
tion becomes better and, once again, close agreement between the two approaches is found, 
supporting our argument put forward in the derivation of the kinetic approach in section II. 
At later times differences become apparent, which may be attributed to the fact that the 
ion relaxation is considerably disturbed by recombination and ionization events leading to 
sudden local changes of the charge density, which is not taken into account by the kinetic 
model. 

Furthermore, according to both approaches, the correlation energy and the thermal ki- 
netic energy of the ions are almost identical roughly to the time where both curves reach 
their maximum values, showing that the total correlation energy is completely converted 
into thermal kinetic energy of the ions, as expressed by Eq. (j!5cj) . At later times, this addi- 
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FIG. 9: Correlation energy (solid line: H-MD simulation, dot-dashed: kinetic result) and thermal 
ion kinetic energy (dashed: H-MD, dotted: kinetic). Initial-state parameters are the same as in 
Fig. ^ The inset shows the ion thermal energy in the early stage of the relaxation process with its 
characteristic transiently oscillatory behavior. 

tional kinetic energy is transferred to the outward directed motion of the ions, leading to an 
indirect enhancement of the plasma expansion by the development of IC and to adiabatic 
cooling of the ions. Therefore, the thermal ion kinetic energy starts to deviate from the 
correlation energy as the plasma expansion sets in. 



IV. CONCLUSIONS 



In summary, we have presented two different theoretical approaches for the simulation of 
ultracold neutral plasmas. First, we have introduced a simple kinetic model along the lines of 
[lol | , and we have shown how to include a description of IC into the model in an approximate 
way. Moreover, we have developed a hybrid molecular dynamics approach which allows for 
an accurate description of the strongly coupled ion motion on microsecond timescales, by 
treating the electronic component as a fluid using an adiabatic approximation while the ions 
are fully accounted for on an MD level. 

n 

Supporting our results from j!3lj . both methods show that the inclusion of IC enhances 
the number of recombined Rydberg atoms by a few percent, but only slightly affects the 
macroscopic expansion dynamics of the plasma itself. As we have shown, this is due to the 
fact that mainly the population of very highly excited states is increased if IC are taken 
into account, which have a small binding energy and therefore hardly influence the electron 
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temperature. 

By comparison of the two methods, we could show that the simple kinetic description 
adequately describes the evolution of global, i.e. spatially averaged, plasma observables. 
Thus, the kinetic model, which allows for a much faster computation, may be used to 
quickly and efficiently scan the vast space of initial-state parameters, e.g. in order to obtain 
a "phase diagram" for Rydberg gas / plasma systems |36j|. Moreover, it permits extending 
the description of ultracold plasmas to a parameter range where the plasma is so large 
that the number of particles (iVj > 10 6 ) prohibits an MD simulation. Maybe even more 
importantly, the simple kinetic equations give additional insight into the dynamics beyond 
that possible on the basis of MD simulations by providing simple evolution equations for 
the macroscopic parameters describing the plasma state. 

On the other hand, spatially resolved quantities such as ionic density, ion velocities or 
local temperature show deviations from the behavior predicted by the simple kinetic model. 
However, the developed H-MD approach provides a powerful method for the study of these 
quantities, and for the detailed description of the relaxation dynamics of the strongly coupled 
ions on a microsecond timescale. Moreover, it permits the st udy of scenarios where the ions 
are so strongly coupled that Coulomb crystallization occurs |l6| , which cannot be described 
by the kinetic model. 
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APPENDIX A: DERIVATION OF THE CORRELATION FORCE 

In this section, the approximation Eq. (|12J) for Fy is derived. We start from Eq. (J7j) 



where the explicit expression tp^ = e 2 /|r — r'| for the inter-ionic Coulomb potential has been 
inserted. In general, the correlation function g$ is a function of both coordinates r and r'. 




(Al) 
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However, in the case of a homogeneous density, ga depends only on the interparticle distance 
x = |r — r'|. Since the relevant property which distinguishes the two points r and r' is the 
corresponding density (from the way the plasma is created, no other differences, e.g. that 
part of the plasma would be in a state with equilibrium correlations while a different part 
would be totally uncorrelated, are apparent), it seems a reasonable approximation to assume 
that the space dependence of the correlation function enters only via the densities at the 
respective coordinates [37]. Hence, we write the correlation function as 

^(r,rOwffii|>i(r),Pi(r , ),|r-r , |] . (A2) 

With the substitution r' = r + x, Eq. (|A1|) becomes 

F a = -e 2 J pi(r + x)cfa [#(r), #(r + x), x] rfx . (A3) 

Since the correlation function rapidly decreases for distances x larger than the correlation 
length A c , we may restrict the integration in Eq. (|A3|) to a sphere with a radius of approx- 
imately A c . If the plasma density does not vary strongly on the scale of the correlation 
length, we may use a linear Taylor expansion of the density 

Pi (r + x) « pi (r) + x • Vpi (r) (A4) 

and the correlation function 



/ / \ / x , dgn(pi, p[, x) 

g n (p u p b x) = gii{p h p h x) + - 



dp[ 



(x-VpO , (A5) 



p[=p 



where p\ = pi(r) and p[ = pi(r + x). 

Substitution of Eq. (|A4|) and Eq. (jA5|) into Eq. ()A3|) and keeping only terms up to linear 
order in x ■ Vp; yields 



x 



F a = -e / -r pi ga(pi,x) rfx + / — ga(p h x) (x • Vp ; ) rfx+- / — p 4 



1 f x dgn(pi,x) 
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dpi 
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(A6) 



where ga(pi,x) = g(pu p[,x)\ p , =p . and we have used the relation 
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which follows from the symmetry of ga under particle exchange, i.e., under exchange of p\ 
and p[. Since the integrand of the first integral in Eq. (jA6|) is an odd function in x the first 
term vanishes. The second term yields after some manipulations 



X 6 

— ga(p h x) (x ■ Vpi) dx = — Vpi 



dx . 



x 



Analogously, the third term can be written as 



- 1 x~z pi 



dga(pi,x) 
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6 dpi 



d f 9n(pi,x) 
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which together leads to 



Fh = -- Vpi 
o 



dpi 
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Pi 
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x 



dx 
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(A8) 



(A9) 



(A10) 



Finally, substitution of the definition of the correlation energy u as given by Eq. (jlHj) yields 



the result Eq. (JT^J). An analogous calcu 
Eq. leads to the familiar LDA result 



ation for the expression of the correlation energy 
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